* Independently-generated stationary outcome and unitroot predictor DGP

* Set working directory
cd "/Users/alikagalwala/Library/CloudStorage/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/Stationary_UnitRoot"

* Run the file with the simulation program
do "Stationary_UnitRoot-DGP.do"

// Monte Carlo Analyses
* Create a matrix to store results
matrix mc_summary_stats = J(35, 73, .)
matrix colnames mc_summary_stats = t phi_y phi_x /// mc identifiers 
b1_stat b1_stat_t1 /// * static model
b1_diff b1_diff_t1 /// * differences model
phi_pa phi_pa_t1 b1_pa b1_pa_t1 lrm_pa lrm_pa_t1  /// * partial adjustment model
phi_ds phi_ds_t1 b2_ds b2_ds_t1  /// * dead start model
phi_gecm phi_gecm_t1 b1_gecm b1_gecm_t1 b2_gecm b2_gecm_t1 lrm_gecm lrm_gecm_t1  /// * GECM
phi_adl phi_adl_t1 b1_adl b1_adl_t1 b2_adl b2_adl_t1 lrm_adl lrm_adl_t1  /// *ADL
phi1_adl22 phi1_adl22_t1 phi2_adl22 phi2_adl22_t1 b1_adl22 b1_adl22_t1 b2_adl22 b2_adl22_t1 b3_adl22 b3_adl22_t1 lrm_adl22 lrm_adl22_t1  /// * ADL22
phi1_adl33 phi1_adl33_t1 phi2_adl33 phi2_adl33_t1 phi3_adl33 phi3_adl33_t1 b1_adl33 b1_adl33_t1 b2_adl33 b2_adl33_t1 b3_adl33 b3_adl33_t1 b4_adl33 b4_adl33_t1 lrm_adl33 lrm_adl33_t1 /// *ADL33
ftest_y_adl33_t1 ftest_y_adl33_power ftest_x_adl33_t1 ///
ftest_y_adl22_t1 ftest_y_adl22_power ftest_x_adl22_t1 phi_pa_power phi_ds_power ///
phi_gecm_power phi_adl_power phi1_adl22_power phi1_adl33_power

* Create a counter to store results in the matrix
local j = 1

foreach t in 25 50 75 100 150 500 1000 {
	foreach phi_y in 0 0.2 0.4 0.6 0.8 {
		foreach phi_x in 1 {
					parallel sim, expr(b1_stat = r(b1_stat) b1_stat_t1 = r(b1_stat_t1) ///
					b1_diff = r(b1_diff) b1_diff_t1 = r(b1_diff_t1) ///
					phi_pa = r(phi_pa) phi_pa_t1 = r(phi_pa_t1) b1_pa = r(b1_pa) b1_pa_t1 = r(b1_pa_t1) lrm_pa = r(lrm_pa) lrm_pa_t1 = r(lrm_pa_t1) ///
					phi_ds = r(phi_ds) phi_ds_t1 = r(phi_ds_t1) b2_ds = r(b2_ds) b2_ds_t1 = r(b2_ds_t1) ///
					phi_gecm = r(phi_gecm) phi_gecm_t1 = r(phi_gecm_t1) b1_gecm = r(b1_gecm) b1_gecm_t1 = r(b1_gecm_t1) b2_gecm = r(b2_gecm) b2_gecm_t1 = r(b2_gecm_t1) lrm_gecm = r(lrm_gecm) lrm_gecm_t1 = r(lrm_gecm_t1) ///
					phi_adl = r(phi_adl) phi_adl_t1 = r(phi_adl_t1) b1_adl = r(b1_adl) b1_adl_t1 = r(b1_adl_t1) b2_adl = r(b2_adl) b2_adl_t1 = r(b2_adl_t1) lrm_adl = r(lrm_adl) lrm_adl_t1 = r(lrm_adl_t1) ///
					phi1_adl22 = r(phi1_adl22) phi1_adl22_t1 = r(phi1_adl22_t1) phi2_adl22 = r(phi2_adl22) phi2_adl22_t1 = r(phi2_adl22_t1) b1_adl22 = r(b1_adl22) b1_adl22_t1 = r(b1_adl22_t1) b2_adl22 = r(b2_adl22) b2_adl22_t1 = r(b2_adl22_t1) b3_adl22 = r(b3_adl22) b3_adl22_t1 = r(b3_adl22_t1) lrm_adl22 = r(lrm_adl22) lrm_adl22_t1 = r(lrm_adl22_t1) ///
					phi1_adl33 = r(phi1_adl33) phi1_adl33_t1 = r(phi1_adl33_t1) phi2_adl33 = r(phi2_adl33)  phi2_adl33_t1 = r(phi2_adl33_t1)phi3_adl33 = r(phi3_adl33) phi3_adl33_t1 = r(phi3_adl33_t1) b1_adl33 = r(b1_adl33) b1_adl33_t1 = r(b1_adl33_t1) b2_adl33 = r(b2_adl33) b2_adl33_t1 = r(b2_adl33_t1) b3_adl33 = r(b3_adl33) b3_adl33_t1 = r(b3_adl33_t1) b4_adl33 = r(b4_adl33) b4_adl33_t1 = r(b4_adl33_t1) lrm_adl33 = r(lrm_adl33) lrm_adl33_t1 = r(lrm_adl33_t1) ftest_y_adl33_t1 = r(ftest_y_adl33_t1) ftest_y_adl33_power = r(ftest_y_adl33_power) ///
					ftest_x_adl33_t1 = r(ftest_x_adl33_t1) ftest_y_adl22_t1 = r(ftest_y_adl22_t1) ftest_y_adl22_power = r(ftest_y_adl22_power) ftest_x_adl22_t1 = r(ftest_x_adl22_t1) phi_pa_power = r(phi_pa_power) phi_ds_power = r(phi_ds_power) ///
					phi_gecm_power = r(phi_gecm_power) phi_adl_power = r(phi_adl_power) ///
					phi1_adl22_power = r(phi1_adl22_power) phi1_adl33_power = r(phi1_adl33_power)) ///
					reps(1000) noisily seeds(12344 12344 12344 12344 12344 12344 12344 12344): NG, obs(`t') phi_y(`phi_y') phi_x(`phi_x')
			
					matrix mc_summary_stats[`j', 1] = `t'
					matrix mc_summary_stats[`j', 2] = `phi_y'
					matrix mc_summary_stats[`j', 3] = `phi_x'
			
					* Static model
					qui su b1_stat
					matrix mc_summary_stats[`j', 4] = r(mean)
					qui su b1_stat_t1
					matrix mc_summary_stats[`j', 5] = r(mean)
					
					* Differences model
					qui su b1_diff
					matrix mc_summary_stats[`j', 6] = r(mean)
					qui su b1_diff_t1
					matrix mc_summary_stats[`j', 7] = r(mean)

					* PA model
					qui su phi_pa
					matrix mc_summary_stats[`j', 8] = r(mean)
					qui su phi_pa_t1
					matrix mc_summary_stats[`j', 9] = r(mean)
					qui su b1_pa
					matrix mc_summary_stats[`j', 10] = r(mean)
					qui su b1_pa_t1
					matrix mc_summary_stats[`j', 11] = r(mean)
					qui su lrm_pa
					matrix mc_summary_stats[`j', 12] = r(mean)
					qui su lrm_pa_t1
					matrix mc_summary_stats[`j', 13] = r(mean)
					
					* DS model
					qui su phi_ds
					matrix mc_summary_stats[`j', 14] = r(mean)
					qui su phi_ds_t1
					matrix mc_summary_stats[`j', 15] = r(mean)
					qui su b2_ds
					matrix mc_summary_stats[`j', 16] = r(mean)
					qui su b2_ds_t1
					matrix mc_summary_stats[`j', 17] = r(mean)
					
					* GECM model
					qui su phi_gecm
					matrix mc_summary_stats[`j', 18] = r(mean)
					qui su phi_gecm_t1
					matrix mc_summary_stats[`j', 19] = r(mean)
					qui su b1_gecm
					matrix mc_summary_stats[`j', 20] = r(mean)
					qui su b1_gecm_t1
					matrix mc_summary_stats[`j', 21] = r(mean)
					qui su b2_gecm
					matrix mc_summary_stats[`j', 22] = r(mean)
					qui su b2_gecm_t1
					matrix mc_summary_stats[`j', 23] = r(mean)
					qui su lrm_gecm
					matrix mc_summary_stats[`j', 24] = r(mean)
					qui su lrm_gecm_t1
					matrix mc_summary_stats[`j', 25] = r(mean)
			
					* ADL model
					qui su phi_adl
					matrix mc_summary_stats[`j', 26] = r(mean)
					qui su phi_adl_t1
					matrix mc_summary_stats[`j', 27] = r(mean)
					qui su b1_adl
					matrix mc_summary_stats[`j', 28] = r(mean)
					qui su b1_adl_t1
					matrix mc_summary_stats[`j', 29] = r(mean)
					qui su b2_adl
					matrix mc_summary_stats[`j', 30] = r(mean)
					qui su b2_adl_t1
					matrix mc_summary_stats[`j', 31] = r(mean)
					qui su lrm_adl
					matrix mc_summary_stats[`j', 32] = r(mean)
					qui su lrm_adl_t1
					matrix mc_summary_stats[`j', 33] = r(mean)
			
					
					* ADL (2,2)
					qui su phi1_adl22
					matrix mc_summary_stats[`j', 34] = r(mean)
					qui su phi1_adl22_t1
					matrix mc_summary_stats[`j', 35] = r(mean)
					qui su phi2_adl22
					matrix mc_summary_stats[`j', 36] = r(mean)
					qui su phi2_adl22_t1
					matrix mc_summary_stats[`j', 37] = r(mean)
					qui su b1_adl22
					matrix mc_summary_stats[`j', 38] = r(mean)
					qui su b1_adl22_t1
					matrix mc_summary_stats[`j', 39] = r(mean)
					qui su b2_adl22
					matrix mc_summary_stats[`j', 40] = r(mean)
					qui su b2_adl22_t1
					matrix mc_summary_stats[`j', 41] = r(mean)
					qui su b3_adl22
					matrix mc_summary_stats[`j', 42] = r(mean)
					qui su b3_adl22_t1
					matrix mc_summary_stats[`j', 43] = r(mean)
					qui su lrm_adl22
					matrix mc_summary_stats[`j', 44] = r(mean)
					qui su lrm_adl22_t1
					matrix mc_summary_stats[`j', 45] = r(mean)
			
					
					* ADL(3,3) model
					qui su phi1_adl33
					matrix mc_summary_stats[`j', 46] = r(mean)
					qui su phi1_adl33_t1
					matrix mc_summary_stats[`j', 47] = r(mean)
					qui su phi2_adl33
					matrix mc_summary_stats[`j', 48] = r(mean)
					qui su phi2_adl33_t1
					matrix mc_summary_stats[`j', 49] = r(mean)
					qui su phi3_adl33
					matrix mc_summary_stats[`j', 50] = r(mean)
					qui su phi3_adl33_t1
					matrix mc_summary_stats[`j', 51] = r(mean)
					qui su b1_adl33
					matrix mc_summary_stats[`j', 52] = r(mean)
					qui su b1_adl33_t1
					matrix mc_summary_stats[`j', 53] = r(mean)
					qui su b2_adl33
					matrix mc_summary_stats[`j', 54] = r(mean)
					qui su b2_adl33_t1
					matrix mc_summary_stats[`j', 55] = r(mean)
					qui su b3_adl33
					matrix mc_summary_stats[`j', 56] = r(mean)
					qui su b3_adl33_t1
					matrix mc_summary_stats[`j', 57] = r(mean)
					qui su b4_adl33
					matrix mc_summary_stats[`j', 58] = r(mean)
					qui su b4_adl33_t1
					matrix mc_summary_stats[`j', 59] = r(mean)
					qui su lrm_adl33
					matrix mc_summary_stats[`j', 60] = r(mean)
					qui su lrm_adl33_t1
					matrix mc_summary_stats[`j', 61] = r(mean)  
					
					* F-tests
					su ftest_y_adl33_t1
					matrix mc_summary_stats[`j', 62] = r(mean) 
					su ftest_y_adl33_power
					matrix mc_summary_stats[`j', 63] = r(mean) 
					su ftest_x_adl33_t1
					matrix mc_summary_stats[`j', 64] = r(mean) 
					su ftest_y_adl22_t1
					matrix mc_summary_stats[`j', 65] = r(mean) 
					su ftest_y_adl22_power
					matrix mc_summary_stats[`j', 66] = r(mean) 
					su ftest_x_adl22_t1
					matrix mc_summary_stats[`j', 67] = r(mean) 
					
					* Power
					su phi_pa_power
					matrix mc_summary_stats[`j', 68] = r(mean)
					su phi_ds_power
					matrix mc_summary_stats[`j', 69] = r(mean)
					su phi_gecm_power
					matrix mc_summary_stats[`j', 70] = r(mean)
					su phi_adl_power
					matrix mc_summary_stats[`j', 71] = r(mean)
					su phi1_adl22_power
					matrix mc_summary_stats[`j', 72] = r(mean)
					su phi1_adl33_power
					matrix mc_summary_stats[`j', 73] = r(mean)
					
					* Iterate the counter
					local j = `j' + 1
			}
		}
}
	

// Store MC matrix results as a dataframe
clear
set obs 35
svmat double mc_summary_stats, names(col)
save stationary_unitroot.dta, replace

* Bias calculations
gen b1_bias_stat = b1_stat - 0
gen phi_bias_pa = phi_pa - phi_y
gen b1_bias_pa = b1_pa - 0
gen phi_bias_ds = phi_ds - phi_y
gen b2_bias_ds = b2_ds - 0
gen phi_bias_adl = phi_adl - phi_y
gen b1_bias_adl = b1_adl - 0
gen b2_bias_adl = b2_adl - 0
gen phi_bias_gecm = phi_gecm - (phi_y - 1)
gen b1_bias_gecm = b1_gecm - 0
gen b2_bias_gecm = b2_gecm - 0
gen phi1_bias_adl22 = phi1_adl22 - phi_y
gen b1_bias_adl22 = b1_adl22 - 0
gen b2_bias_adl22 = b2_adl22 - 0
gen b1_bias_diff = b1_diff - 0
gen phi1_bias_adl33 = phi1_adl33 - phi_y
gen b1_bias_adl33 = b1_adl33 - 0
gen b2_bias_adl33 = b2_adl33 - 0
gen lrm_bias_pa = lrm_pa - 0
gen lrm_bias_adl = lrm_adl - 0
gen lrm_bias_gecm = lrm_gecm - 0
gen lrm_bias_adl22 = lrm_adl22 - 0
gen lrm_bias_adl33 = lrm_adl33 - 0
gen phi2_bias_adl22 = phi2_adl22 - 0
gen phi2_bias_adl33 = phi2_adl33 - 0
gen phi3_bias_adl33 = phi3_adl33 - 0
gen b3_bias_adl22 = b3_adl22 - 0
gen b3_bias_adl33 = b3_adl33 - 0
gen b4_bias_adl33 = b4_adl33 - 0

save stationary_unitroot.dta, replace
